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Abstract 

In this research paper the analytical expression for the 
concentrations of the reaction-diffusion equations in steady 
state Brusselator model is discussed. The approximate 
analytical solution of the steady state non-linear boundary 
value problem is derived by using the Homotopy 
perturbation method. The Homotopy-perturbation method 
has the advantage of being more concise for analytical and 
numerical functions. This method can be easily extended to 
find the solution of the other strongly non-linear initial and 
boundary value problems in engineering and chemical 
sciences. 
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Introduction 

Systems of ordinary differential equations (ODEs) 
arise in mathematical models throughout science and 
engineering. When an explicit condition (or 

conditions) that a solution must satisfy is specified at 
one value of the independent variable, usually its 
lower bound, this is referred to as an initial value 
problem (I VP). When the conditions to be satisfied 
occur at more than one value of the independent 
variable, this is referred to as a boundary value 
problem (BVP). If there are two values of the 
independent variable at which conditions are 
specified, then this is a two-point boundary value 
problem (TPBVP). TPBVPs occur in a wide variety of 
problems, including the modeling of chemical 
reactions, heat transfer, and diffusion. They are also of 
interest in optimal control problems. For a typical 
TPBVP, when expressed as system of first-order ODEs, 
there are one or more unknown initial states. Solving 
the TPBVP is equivalent to finding these initial states. 
In other cases, the initial states are known, but there 
are unknown model parameters, and solving the 
TPBVP requires finding values for these parameters. 


There may also be a combination of unknown initial 
states and parameters that must be determined to 
solve theTPBVP. A TPBVP may have no solution, a 
single solution, or multiple solutions. There are many 
techniques available for the numerical solution of 
TPBVPs for ODEs [Heath (2002), Chen et. al (2005 & 
2006)]. 

Reaction-diffusion (RD) systems arise frequently in the 
study of chemical and biological phenomena and are 
naturally modeled by partial differential equations. 
Finding accurate and efficient methods for solving 
nonlinear system of PDEs has long been an active 
research undertaking. In this analysis, we study the 
reaction-diffusion Brusselator model which models a 
chemical reaction diffusion process by homotopy 
perturbation method. Wazwaz (2000) used the 
Adomian decomposition method to handle the 
reaction-diffusion Brusselator model. In recent years, 
much attention has been given to the study of the 
Homotopy-perturbation method (HPM) for solving a 
wide range of problems whose mathematical models 
yield differential equation or system of differential 
equations. 

Mathematical Formulations of the Problem 

The steady state Brusselator model is a reaction- 
diffusion system in fixed, bounded and one 
dimensional medium. The reaction scheme of the 
system is [Lin et. al (2008)] 

A <-> V 

2X+Y <->3X 

B + X^Y + D 

X <-> E 

The system is open, allowing the concentrations of the 
reactants and to be kept at a constant level 
throughout the reaction/dif fusion medium. In the 
event of all forward rate constants set to one and all 
reverse rate constants set to zero, the steady-state 
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concentration profiles for X and Y are described are 
by 

p 2 v 

D x — r = (B + l)X-A-X 2 Y (1) 

d 2 Y , 

D Y — — = X 2 Y - B X (2) 

dz 2 

Where A , B , X and Y are real and denote the 
concentrations of the corresponding chemical species, 
D x and D Y are the diffusion coefficients for X and 
Y respectively. The spatial coordinates in the 
reaction/diffusion medium is the independent 
variable ze[0,L] . The corresponding boundary 

conditions are 

X(0) = X(L) = A and 7(0) = 7(L) = B/A (3) 

Solution of the Boundary Value Problem 
Using the Homotopy Perturbation Method 

Linear and non-linear phenomena are of fundamental 
importance in various fields of science and 
engineering. Most models of real - life problems are 
still very difficult to solve. Therefore, approximate 
analytical solutions such as Homotopy perturbation 
method (HPM) [Ghori et. al (2007), Ozis et. al(2007), Li 
et. al(200 7), Mousa et. al(2008). He (1999 and 2003), 
Ariel (2010), Ananthaswamy et. al(2012) and Subha et. 
al(2014) ] were introduced. This method is the most 
effective and convenient ones for both linear and non- 
linear equations. Perturbation method is based on 
assuming a small parameter. The majority of non- 
linear problems, especially those having strong non- 
linearity, have no small parameters at all and the 
approximate solutions obtained by the perturbation 
methods, in most cases, are valid only for small values 
of the small parameter. Generally, the perturbation 
solutions are uniformly valid as long as a scientific 
system parameter is small. However, we cannot rely 
fully on the approximations, because there is no 
criterion on which the small parameter should exist. 
Thus, it is essential to check the validity of the 
approximations numerically and/or experimentally. 
To overcome these difficulties, HPM have been 
proposed recently. 

Recently, many authors have applied the Homotopy 
perturbation method (HPM) to solve the nonlinear 
boundary value problems in physics and engineering 
sciences [Ghori et. al (2007), Ozis et. al(2007), Li et. 


al(2007), Mousa et. al(2008)]. Recently this method is 
also used to solve some of the non-linear problem in 
physical sciences [He (1999 and 2003)]. This method is 
a combination of Homotopy in topology and classic 
perturbation techniques. Ji-Huan He used to solve the 
Lighthill equation (1999), the Diffusion equation (2003) 
and the Blasius equation [He (2003) and Ariel (2010)]. 
The HPM is unique in its applicability, accuracy and 
efficiency. The HPM uses the imbedding parameter p 
as a small parameter, and only a few iterations are 
needed to search for an asymptotic solution. 

Appendix A: Basic Concept of the 

Homotopy Perturbation Method [He (1999 
and 2003)] 

To explain this method, let us consider the following 
function: 

D 0 (u) — f (r) = 0, reQ (A.l) 

with the boundary conditions of 

Pu 

B 0 (u, — ) = 0, reT (A.2) 

on 

where D 0 is a general differential operator, B 0 is a 
boundary operator, /(r) is a known analytical 
function and T is the boundary of the domain Q . In 
general, the operator D 0 can be divided into a linear 
part L and a non-linear part N . The eqn.(A.l) can 
therefore be written as 

L(u) + N(u)-f(r) = 0 (A.3) 

By the Homotopy technique, we construct a 
Homotopy v(r, p ) : Q x [0, 1] — » 91 that satisfies 

H(y,p) = (l-p)[L(v)-L(u 0 )] + p[D o (v)-f(r)]=0 (A.4) 
H(v , p) = L(y) - L(u 0 ) + pL(u 0 ) + p[N(v) - f(r )] = 0 (A.5) 
where p e [0, 1] is an embedding parameter, and u 0 is 
an initial approximation of the eqn.(A.l) that satisfies 
the boundary conditions. From the eqns. (A.4) and 
(A.5), we have 

H(y , 0) = L(v) - L(u 0 ) = 0 (A.6) 

H (v, Y) = D 0 (v) — f (r) = 0 (A.7) 

When p= 0, the eqns. (A.4) and (A.5) become linear 
equations. When p =1, they become non-linear 
equations. The process of changing p from zero to 
unity is that of L(v)-L(u 0 ) = 0to D 0 (v)-f(r) =0 . We 
first use the embedding parameter p as a small 
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parameter and assume that the solutions of the eqns. 
(A.4) and (A.5) can be written as a power series in p : 

v = v 0 + pv x + p 2 v 2 + ... (A. 8) 

Setting p = 1 results in the approximate solution of 
the eqn. (A.l): 


u = lim v = v 0 + Vj + v 2 + ... 
p->i 


(A.9) 


This is the basic idea of the HPM. 

Appendix B: Solution of the Non-linear Eqns. 
(1)-(3) by Using the Homotopy 
Perturbation Method 


In this Appendix, we indicate how the eqns. (4) and (5) 
is derived in this paper. To find the solution of the 
eqns.(l)-(3) ,we construct the Homotopy as follows: 


C l ~p ) 


d 2 X (B + ])X 


+ P 


dz 2 

<px_ 

dZ 2 


Dy 


\B + \)X-A-X 2 Y ^ 


D v 


=0 


(B.l) 


(1 ~P) 


j 2 f 

dz 2 


+ p 


d 2 Y ( X 2 Y-BX ) 


dz 2 


D r 


=0 


(B.2) 


The analytical solution of the eqns. (B.l) and (B.2) are 


X = X 0 +pX l +p 2 X 2 +. 
Y = Y 0 + pY l + p 2 Y 2 +.... 


Substituting the eqns.(B.3) and (B.4) 
eqns. (B.l) and (B.2) respectively we get 


(B-3) 
(B.4) 
into the 


0-/0 


i 2 (x 0 + py + p 2 x 2 +....) 

dz 2 

(g + l)(X 0 + pX 1+ p 2 X 2 +....) 

d y 


+p 


t/ 2 (X 0 +pX 1 + p 2 X 2 +...) 

dz 2 

f (B + 1)(X 0 + pX ] + p 2 X 2 +....) - A 
(X 0 -pX 1 + p 2 X 2 +....) 2 


(F 0 + pF 1 + p 2 F 2 +...) 

Dy 


(B.5) 


=0 


0-/0 


+ P 


f(Y^ + pY l + PY 2 + ) 

dz 2 

+ PY\ + P^2 + ) 

dz 2 

\x 0 + p x 1 + p 2 x 2 +....) 2 

<Y () + pY y + p 2 Y 2 + ) 

-g(X 0 + pX 1+; p 2 X 2 +....) 
D v 


(B.6) 


= 0 


V y_ 

Comparing the coefficients of like powers of p in the 
eqns. (B.5) and (B.6) we get 

72, 


o . d z X Q 

P ' * 2 


dz z 


^ {B + Y)XP 

v D x J 


= 0 


o.A =0 


P ' dz 2 


, d 2 X, 

7 : — ? 

dz 2 


(B + 1)XP 


a+x 2 f 0 ) 

K D X ) 

1 

l D x J 


dz 2 


BX^-XljP 

Dy 


= 0 


(B.7) 

(B.8) 
= 0 (B.9) 

(B.10) 


The initial approximations are as follows: 

X 0 (0) = X 0 (L) = A and F 0 (0) = Y 0 (L) = B / A (B.ll) 

X f (0) = X t (L) = 0 and 
Y 0 (0) = Y 0 (L) = B/ A,i = 1,2,3 


(B.12) 


Solving the eqns.(B.7)-(B.10) and using the boundary 
condition eqn. (B.ll) and (B.12), we obtain the 
following results: 




2sinh(&L) 


-(l-e kL )e~ kz ] (B.13) 


Y-* 

0_ A 


(B.14) 


v = 


A V(1 - e - kL )e kz -(l-e kL ) e - kz ]- 

2 sinh(£L) J L J 


1+ 

f A ) 

D x y 

J 


AB 


(1 -e~ kL ) 2 e 2kz +{\-e KL ) z e 
+6(\-e kL ) (l-e~ kL ) 


kL \2 -2k z 


12 k 2 D x sinh 2 (kL) 


(B.15) 
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' AB ' 


' AB((l-e kL ) 2 +(l-e~ kL ) 2 y 

[Dyk 2 j 


16 D y k 2 sinh 2 (k L) 

v y 


AB{l-e kL ){l-e~ kL ){Lz-z 2 )\ 


V 

A D Y sinh 2 {kL) y 



(AB(i 

[l-e~ kL )e kz -{l-e kL )e~ ki 

: )1 


V 

2 D y sinh 2 (A: L) 

y 


' AB({ 

l-e~ kL ) 2 e 2kz +{l-e kL ) 2 e 

-2k z 

V 

16k 2 D y sinh 1 {kL) 


y 


where k and k x are defined by the eqns. (6) and (7) 
respectively. 

According to the HPM, we can conclude that 

X = limX(z) = X 0 + X! (B.17) 

i 

Y = ]imY(z) = Y 0 +Y 1 (B.18) 

p^i 


After putting the eqns. (B.13) and (B.15) into an 
eqn.(B.17), and (B.14) and (B.16) into an eqn. (B.18), we 
obtain the solution in the text eqns. (4) and (5). 


The approximate analytical solution of eqns.(l)-(3) 
using the Homotopy perturbation method is given by 




k x +A A 


2 sinh(^L) 

^ A '' 


({l-e~ kL )e kz -{l-e kL )e~ kz ) 


k 2 D 

vA u x y 


(4) 


( r {l-e~ kL ) 2 e 2kz +{l -e kL fe~ 2k ^ 


AB 


+6 (l-e kL ) (l-e~ kL ) 


12 k 2 D x sinh 2 (kL) 


V 


y 


' AB ' 


f AB({l-e kL ) 2 +{l-e~ kL ) 2 y 

{Dyk 2 J 


16 D Y k 2 sinh 2 (/: L) 

\ y 


AB{l-e kL ){l-e~ kL ){Lz-Z 2 )\ 


V 

' AB[{ 

4 D Y sinh 2 (kL) J 

l-e~ kL )e kz -{l-e kL )e~ kz y 


V 

2 D y sinh 2 {kL) 


' AB({ 

1 - e~ kL ) 2 e 2k z +{l-e kL f e~ 2kz y 

V 

16k 2 D y sinh 2 {k L) 

y 


(5) 


where 


k = 


\B + 1 ) 

D x 


K = 


AB 


v +6(l-/ L )(l- g - tL ) y 


12 k 2 D x sinh 2 (kL) 


(6) 


(7) 


Results and Discussions 

Figure 1 and 2 shows the dimensionless concentration 
X versus the dimensionless distance z. From the Fig. 
1 (a)- (g), it is clear that, when the dimensionless 
constant L increases, the corresponding 
concentration X (z) decreases in some fixed values of 
the dimensionless parameters D x , D Y , A and B. From 
Fig.2 (a)-(b), it is noted that when the diffusion 
coefficient D x increases, the corresponding 
concentrations X(z) also increases in some fixed 
values of the dimensionless parameters and for 
various ranges of the dimensionless constant L . 


(a) 
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(c) 



(d) 



(e) 



(f) 
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(g) 



FIGURE 1: DIMENSIONLESS CONCENTRATION X VERSUS 
THE DIMENSIONLESS DISTANCE z • THE CONCENTRATION 
X(z) WERE COMPUTED USING THE EQN. (4) IN SOME FIXED 
VALUES OF THE DIMENSIONLESS PARAMETERS 
D x , D Y , A and B , WHEN (A) L = 0.05,(B) L = 0.10 / (C) L = 0.15, 
(D) L = 0.20, (E) L = 0.25, (F) L = 0.30 AND (G) L = 0.35 . 


(a) 



(b) 



FIGURE 2: DIMENSIONLESS CONCENTRATION X VERSUS 
THE DIMENSIONLESS DISTANCE z • THE CONCENTRATION 
X( z ) WERE COMPUTED USING THE EQN. (4) IN SOME FIXED 
VALUES OF THE DIMENSIONLESS PARAMETERS D r , A and B 
AND VARIOUS VALUES OF THE DIFFUSION COEFFICIENT D x 
WHEN (A) L = 0.10 AND (B) L = 0.35 . 
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(a) 



Dime monies s distance z 


(b) 



(c) 



(d) 



Advances in Chemical Science Volume 3 Issue 3, November 2014 


(e) 




(g) 



FIGURE 3: DIMENSIONLESS CONCENTRATION Y VERSUS 
THE DIMENSIONLESS DISTANCE z . THE CONCENTRATION 
Y{z) WERE COMPUTED USING THE EQN. (4) IN SOME FIXED 
VALUES OF THE DIMENSIONLESS PARAMETERS 
D x , D Y , A and B , WHEN (A) L = 0.05 , (B) L = 0.10,(C) L = 0.15, 
(D) L = 0.20, (E) L - 0.25 , (F) L = 0.30 AND (G) L = 0.35 . 

Figure 3 and 4 shows the dimensionless concentration 
Y versus the dimensionless distance z- From the Fig. 2 
(a)- (g), it is inferred that when the dimensionless 
constant L increases, the corresponding 
concentration Y(z) also increases in some fixed values 
of the dimensionless parameters D x , D Y , A and B. 
From Fig.4 (a)-(b), it is observed that when the 
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diffusion coefficient D x increases, the corresponding 
concentrations Y(z) also increases in some fixed 
values of the dimensionless parameters and for 
various ranges of the dimensionless constant L . 


(a) 




FIGURE 4: DIMENSIONLESS CONCENTRATION X VERSUS 
THE DIMENSIONLESS DISTANCE z • THE CONCENTRATION 
X(z) WERE COMPUTED USING THE EQN. (4) IN SOME FIXED 
VALUES OF THE DIMENSIONLESS PARAMETERS D y , A and B 
AND VARIOUS VALUES OF THE DIFFUSION COEFFICIENT D x 
WHEN (A) L = 0.10 AND (B) L = 0.35 . 

Conclusion 

This investigation is carried out using the Homotopy 
perturbation method (HPM) to approximate the 
solutions of well-known steady state Brusselator 
reaction-diffusion model. The HPM is very simple and 
straightforward. HPM avoids the difficulties arising in 
finding the Adomian polynomials and transformation 
formulas. Very recently, the HPM has been widely 
used to solve many problems in engineering and 
sciences because of it's reliability and the reduction in 
the size of computations. It is shown that the 


Homotopy-perturbation method is a promising tool 
for linear and non-linear system of ordinary and 
partial differential equations. 
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Appendix C: Nomenclature 


Symbo 

1 

Meaning 

A 

Dimensionless constant 

B 

Dimensionless constant 

X 

Steady state concentration profile for 
the chemical specie A 

Y 

Steady state concentration profile for 
the chemical specie B 

z 

Dimensionless distance 

L 

Dimensionless constant 

D x 

Diffussion coefficient for the 
concentration of the chemical specie A 

Dy 

Diffussion coefficient for the 
concentration of the chemical specie B 
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